%Matlab code to generate figures 4 in the manuscript "Putting the
%significance of spectral peaks on the level: implications for the
%1470-yr peak in Greenland δ18O" by P. Huybers, Feb 2022.  This
%file was run on Matlab release R2021b and requires the signal
%processing toolbox. 

%Various statistics that are reported in Sections 4 and 5 (simulated results)

clear;
load GISPsynth.mat x P A s k time opts Vexp;

%Spectral energy within 1/1470 \pm 0.1/1470 for H0
dummy=P(1).init;
pl=find(s>(1/1.470-0.1/1.470) & s<(1/1.470+0.1/1.470));
val=mean(sum(exp(dummy(:,pl))')./sum(exp(dummy)')); %6.56 percent

%Spectral energy within 1/1470 \pm 0.1/1470 for H1
dummy=A(1).init;
pl=find(s>(1/1.470-0.1/1.470) & s<(1/1.470+0.1/1.470));
val=mean(sum(exp(dummy(:,pl))')./sum(exp(dummy)')); %22.41 percent

%H0(1,0)max false rejection rate
alpha=0.01;
conf=log(gaminv(1-alpha,k,1/k));
confb=log(gaminv(1-alpha*2*k/length(time),k,1/k));
H010max=max(P(1).pre');
mean(H010max>confb); %0.6104

%Average F for H_0(1,0)
F=var(P(1).pre')/Vexp; 
mean(F); %2.1494
prctile(F,[2.5 97.5]); %1.3962    3.2461

%significance levels using the results of the simulations (black vertical line in Fig.~4a)
dummy=prctile(H010max,99); %2.1044
load EX_GISP.mat mag1470;
mag1470(1)<dummy; % yes

%H0(1,1)max false rejection rate
alpha=0.01;
conf=log(gaminv(1-alpha,k,1/k));
confb=log(gaminv(1-alpha*2*k/length(time),k,1/k));
H010max=max(P(2).pre');
mean(H010max>confb); %0.0154

%Average F for H_0(1,1)
F=var(P(2).pre')/Vexp; 
mean(F); %1.0270
prctile(F,[2.5 97.5]); %0.6418    1.6151

%Power of the test
for ct=1:3;
    ho=max(P(ct).pre');
    cv=prctile(ho,99);
    ao=max(A(ct).pre');
    power(ct)=mean(ao>cv); %0.6306    0.4995    0.0800
end;
